1 Import data

languages <- read_csv("voicing-effect/stimuli/languages.csv")
## Parsed with column specification:
## cols(
##   speaker = col_character(),
##   language = col_character()
## )
words <- read_csv("voicing-effect/stimuli/nonce.csv")
## Parsed with column specification:
## cols(
##   item = col_integer(),
##   word = col_character(),
##   ipa = col_character(),
##   c1 = col_character(),
##   c1phonation = col_character(),
##   vowel = col_character(),
##   anteropost = col_character(),
##   height = col_character(),
##   c2 = col_character(),
##   c2phonation = col_character(),
##   c2place = col_character(),
##   language = col_character()
## )
columns <- c(
    "speaker",
    "seconds",
    "rec.date",
    "prompt",
    "label",
    "TT.displacement.sm",
    "TT.velocity",
    "TT.velocity.abs",
    "TD.displacement.sm",
    "TD.velocity",
    "TD.velocity.abs"
)

aaa.files <- list.files(
    path = "./voicing-effect/results/ultrasound",
    pattern = "*-tongue-cart.tsv",
    full.names = TRUE
)

tongues <- read_aaa(
    aaa.files,
    columns, 
    na.rm = TRUE
) %>%
    mutate(word = word(prompt, 2)) %>%
    left_join(y = languages) %>%
    left_join(y = words) %>%
    mutate_if(is.character, as.factor) %>%
    group_by(speaker) %>%
    mutate(
        X.re = rescale(X),
        Y.re = rescale(Y)
    ) %>%
    ungroup() %>%
    mutate(
        vowel.ord = ordered(vowel, levels = c("a", "o", "u")),
        c2place.ord = ordered(c2place, levels = c("coronal", "velar")),
        c2phonation.ord = ordered(c2phonation, levels = c("voiceless", "voiced"))
    )
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   speaker = col_character(),
##   rec.date = col_character(),
##   prompt = col_character(),
##   label = col_character()
## )
## See spec(...) for full column specifications.
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   speaker = col_character(),
##   rec.date = col_character(),
##   prompt = col_character(),
##   label = col_character()
## )
## See spec(...) for full column specifications.
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   speaker = col_character(),
##   rec.date = col_character(),
##   prompt = col_character(),
##   label = col_character()
## )
## See spec(...) for full column specifications.
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   speaker = col_character(),
##   rec.date = col_character(),
##   prompt = col_character(),
##   label = col_character(),
##   X_1 = col_character(),
##   Y_1 = col_character()
## )
## See spec(...) for full column specifications.
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   speaker = col_character(),
##   rec.date = col_character(),
##   prompt = col_character(),
##   label = col_character()
## )
## See spec(...) for full column specifications.
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   speaker = col_character(),
##   rec.date = col_character(),
##   prompt = col_character(),
##   label = col_character()
## )
## See spec(...) for full column specifications.
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   speaker = col_character(),
##   rec.date = col_character(),
##   prompt = col_character(),
##   label = col_character()
## )
## See spec(...) for full column specifications.
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   speaker = col_character(),
##   rec.date = col_character(),
##   prompt = col_character(),
##   label = col_character()
## )
## See spec(...) for full column specifications.
## Joining, by = "speaker"
## Joining, by = c("word", "language")
max <- tongues %>%
    filter(label %in% c("max_TT", "max_TD"), vowel != "u") %>%
    arrange(rec.date, fan.line) %>%
    create_event_start("rec.date")

max_it_12 <- max %>%
    filter(speaker %in% c("it01", "it02"))

max_pl_23 <- max %>%
    filter(speaker %in% c("pl02", "pl03"))

max_pl_34 <- max %>%
    filter(speaker %in% c("pl03", "pl04"))

aaa_files_clos <- list.files(
    path = "./voicing-effect/results/ultrasound",
    pattern = "*-tongue-clos-cart.tsv",
    full.names = TRUE
)

tongues_clos <- read_aaa(
    aaa_files_clos,
    columns, 
    na.rm = TRUE
) %>%
    mutate(word = word(prompt, 2)) %>%
    left_join(y = languages) %>%
    left_join(y = words) %>%
    mutate_if(is.character, as.factor) %>%
    group_by(speaker) %>%
    mutate(
        X.re = rescale(X),
        Y.re = rescale(Y)
    ) %>%
    ungroup() %>%
    mutate(
        vowel.ord = ordered(vowel, levels = c("a", "o", "u")),
        c2place.ord = ordered(c2place, levels = c("coronal", "velar")),
        c2phonation.ord = ordered(c2phonation, levels = c("voiceless", "voiced"))
    ) %>%
    arrange(rec.date, fan.line) %>%
    create_event_start("rec.date")
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   speaker = col_character(),
##   rec.date = col_character(),
##   prompt = col_character(),
##   label = col_character(),
##   X_2 = col_character(),
##   Y_2 = col_character(),
##   X_3 = col_character(),
##   Y_3 = col_character(),
##   X_4 = col_character(),
##   Y_4 = col_character(),
##   X_5 = col_character(),
##   Y_5 = col_character(),
##   X_6 = col_character(),
##   Y_6 = col_character()
## )
## See spec(...) for full column specifications.
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   speaker = col_character(),
##   rec.date = col_character(),
##   prompt = col_character(),
##   label = col_character(),
##   X_3 = col_character(),
##   Y_3 = col_character(),
##   X_4 = col_character(),
##   Y_4 = col_character(),
##   X_5 = col_character(),
##   Y_5 = col_character()
## )
## See spec(...) for full column specifications.
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   speaker = col_character(),
##   rec.date = col_character(),
##   prompt = col_character(),
##   label = col_character()
## )
## See spec(...) for full column specifications.
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   speaker = col_character(),
##   rec.date = col_character(),
##   prompt = col_character(),
##   label = col_character()
## )
## See spec(...) for full column specifications.
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   speaker = col_character(),
##   rec.date = col_character(),
##   prompt = col_character(),
##   label = col_character()
## )
## See spec(...) for full column specifications.
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   speaker = col_character(),
##   rec.date = col_character(),
##   prompt = col_character(),
##   label = col_character()
## )
## See spec(...) for full column specifications.
## Joining, by = "speaker"
## Joining, by = c("word", "language")

2 Exploration

max %>%
    filter(c2place == "coronal", language == "italian") %>%
    plot_tongue(geom = "point", alpha = 0.5) +
    aes(colour = c2phonation) +
    facet_grid(vowel ~ speaker)

max %>%
    filter(c2place == "velar", language == "italian") %>%
    plot_tongue(geom = "point", alpha = 0.5) +
    aes(colour = c2phonation) +
    facet_grid(vowel ~ speaker)

max %>%
    filter(c2place == "coronal", language == "polish") %>%
    ggplot(aes(X, Y, colour = c2phonation)) +
    geom_point(alpha = 0.5) +
    facet_grid(vowel ~ speaker) +
    coord_fixed(1.5)

max %>%
    filter(c2place == "coronal", language == "polish") %>%
    ggplot(aes(X.re, Y.re, colour = c2phonation)) +
    geom_point(alpha = 0.5) +
    facet_grid(speaker ~ vowel) +
    coord_fixed(0.7)

max %>%
    filter(c2place == "velar", language == "polish") %>%
    ggplot(aes(X, Y, colour = c2phonation)) +
    geom_point(alpha = 0.5) +
    facet_grid(vowel ~ speaker) +
    coord_fixed(0.7)

max %>%
    filter(c2place == "velar", language == "polish") %>%
    ggplot(aes(X.re, Y.re, colour = c2phonation)) +
    geom_point(alpha = 0.5) +
    facet_grid(speaker ~ vowel) +
    coord_fixed(0.7)

3 Testing

Comparing fixed effects models works best with ML, otherwise you can use fREML with AIC for fixed effects, but if there is an AR1 model cannot use AIC, so need ML.

3.1 Polish

pl.m1 <- bam(
    Y.re ~
        X.re +
        s(X.re, bs = "cr") +
        s(X.re, by = c2phonation.ord, bs = "cr") +
        s(X.re, by = c2place.ord, bs = "cr") +
        s(X.re, by = vowel.ord, bs = "cr") +
        s(X.re, rec.date, bs = "fs", xt = "cr", m = 1, k = 5) +
        s(X.re, speaker, bs = "fs", xt = "cr", m = 1, k = 5),
    data = max_pl_23,
    method = "fREML"
)

rho <- start_value_rho(pl.m1)

pl.m1.ar <- bam(
    Y.re ~
        X.re +
        s(X.re, bs = "cr") +
        s(X.re, by = c2phonation.ord, bs = "cr") +
        s(X.re, by = c2place.ord, bs = "cr") +
        s(X.re, by = vowel.ord, bs = "cr") +
        s(X.re, rec.date, bs = "fs", xt = "cr", m = 1, k = 5) +
        s(X.re, speaker, bs = "fs", xt = "cr", m = 1, k = 5),
    data = max_pl_23,
    method = "fREML",
    rho = rho,
    AR.start = max_pl_23$start.event
)

summary(pl.m1)

pl.m1.null <- bam(
    Y.re ~
        X.re +
        s(X.re, bs = "cr") +
        s(X.re, by = c2place.ord, bs = "cr") +
        s(X.re, by = vowel.ord, bs = "cr") +
        s(X.re, rec.date, bs = "fs", xt = "cr", m = 1, k = 5) +
        s(X.re, speaker, bs = "fs", xt = "cr", m = 1, k = 5),
    data = max_pl_23,
    method = "fREML"
)

compareML(pl.m1, pl.m1.null)
plot_smooth(
    pl.m1,
    view = "X.re",
    plot_all= "c2phonation.ord",
    cond = list("c2place.ord" = "coronal"),
    rug = FALSE
)
plot_diff(
    pl.m1,
    view = "X.re",
    comp = list(c2phonation.ord = c("voiceless", "voiced")),
    cond = list("c2place.ord" = "coronal")
)
plot_gamsd(
    pl.m1,
    view = "X.re",
    comparison = list(c2phonation.ord = c("voiceless", "voiced")),
    conditions = list(c2place.ord = "coronal")
)
plot_gamsd(
    pl.m1.ar,
    view = "X.re",
    comparison = list(c2phonation.ord = c("voiceless", "voiced")),
    conditions = list(c2place.ord = "coronal")
)

3.1.1 PL02

pl02_max <- filter(max, speaker == "pl02", X > -20)

pl02.m1 <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = pl02_max,
    method = "ML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
rho <- start_value_rho(pl02.m1)

pl02.m1.ar <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = pl02_max,
    method = "ML",
    rho = rho,
    AR.start = pl02_max$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
summary(pl02.m1)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ X + s(X, bs = "cr") + s(X, by = c2phonation.ord, bs = "cr") + 
##     s(X, by = c2place.ord, bs = "cr") + s(X, by = vowel.ord, 
##     bs = "cr") + s(X, rec.date, bs = "fs", xt = "cr", m = 1, 
##     k = 5)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.23973    0.59721   5.425 6.45e-08 ***
## X           -0.09449    0.03601  -2.624  0.00876 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                edf  Ref.df       F  p-value    
## s(X)                         7.269   7.773  66.456  < 2e-16 ***
## s(X):c2phonation.ordvoiced   8.003   8.684   6.916 5.04e-07 ***
## s(X):c2place.ordvelar        8.877   8.983 309.316  < 2e-16 ***
## s(X):vowel.ordo              7.803   8.573  21.311  < 2e-16 ***
## s(X,rec.date)              188.143 300.000   9.608  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Rank: 342/343
## R-sq.(adj) =  0.949   Deviance explained = 95.4%
## -ML =   5542  Scale est. = 4.4805    n = 2391
pl02.m1.null <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = pl02_max,
    method = "ML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
compareML(pl02.m1, pl02.m1.null)
## pl02.m1: Y ~ X + s(X, bs = "cr") + s(X, by = c2phonation.ord, bs = "cr") + 
##     s(X, by = c2place.ord, bs = "cr") + s(X, by = vowel.ord, 
##     bs = "cr") + s(X, rec.date, bs = "fs", xt = "cr", m = 1, 
##     k = 5)
## 
## pl02.m1.null: Y ~ X + s(X, bs = "cr") + s(X, by = c2place.ord, bs = "cr") + 
##     s(X, by = vowel.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5)
## 
## Chi-square test of ML scores
## -----
##          Model    Score Edf  Chisq    Df   p.value Sig.
## 1 pl02.m1.null 5557.001  11                            
## 2      pl02.m1 5541.991  13 15.011 2.000 3.027e-07  ***
## 
## AIC difference: -35.85, model pl02.m1 has lower AIC.
plot_gamsd(
    pl02.m1.ar,
    view = "X",
    comparison = list(c2phonation.ord = c("voiceless", "voiced")),
    conditions = list(c2place.ord = "coronal")
)
## Summary:
##  * X : numeric predictor; with 100 values ranging from -19.967800 to 59.219800. 
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * rec.date : factor; set to the value(s): 07/02/2017 16:29:14.

3.1.2 PL03

pl03_max <- filter(max, speaker == "pl03")

pl03.m1 <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = pl03_max,
    method = "fREML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
rho <- start_value_rho(pl03.m1)

pl03.m1.ar <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = pl03_max,
    method = "ML",
    rho = rho,
    AR.start = pl03_max$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
pl03.m1.ar.null <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
#        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = pl03_max,
    method = "ML",
    rho = rho,
    AR.start = pl03_max$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
summary(pl03.m1)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ X + s(X, bs = "cr") + s(X, by = c2phonation.ord, bs = "cr") + 
##     s(X, by = c2place.ord, bs = "cr") + s(X, by = vowel.ord, 
##     bs = "cr") + s(X, rec.date, bs = "fs", xt = "cr", m = 1, 
##     k = 5)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -5.15251    0.84915  -6.068 1.63e-09 ***
## X            0.42071    0.05257   8.004 2.35e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                edf  Ref.df      F  p-value    
## s(X)                         7.433   7.793 52.044  < 2e-16 ***
## s(X):c2phonation.ordvoiced   6.740   7.654  4.817 1.54e-05 ***
## s(X):c2place.ordvelar        8.755   8.938 33.136  < 2e-16 ***
## s(X):vowel.ordo              8.770   8.941 20.139  < 2e-16 ***
## s(X,rec.date)              209.072 235.000 26.599  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Rank: 277/278
## R-sq.(adj) =  0.969   Deviance explained = 97.3%
## fREML = 3995.5  Scale est. = 2.8743    n = 1793
pl03.m1.null <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = pl03_max,
    method = "ML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
compareML(pl03.m1.ar, pl03.m1.ar.null)
## pl03.m1.ar: Y ~ X + s(X, bs = "cr") + s(X, by = c2phonation.ord, bs = "cr") + 
##     s(X, by = c2place.ord, bs = "cr") + s(X, by = vowel.ord, 
##     bs = "cr") + s(X, rec.date, bs = "fs", xt = "cr", m = 1, 
##     k = 5)
## 
## pl03.m1.ar.null: Y ~ X + s(X, bs = "cr") + s(X, by = c2place.ord, bs = "cr") + 
##     s(X, by = vowel.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5)
## 
## Chi-square test of ML scores
## -----
##             Model    Score Edf Chisq    Df p.value Sig.
## 1 pl03.m1.ar.null 3335.866  11                         
## 2      pl03.m1.ar 3330.848  13 5.018 2.000   0.007  **
## Warning in compareML(pl03.m1.ar, pl03.m1.ar.null): AIC is not reliable,
## because an AR1 model is included (rho1 = 0.611727, rho2 = 0.611727).
plot_gamsd(
    pl03.m1.ar,
    view = "X",
    comparison = list(c2phonation.ord = c("voiceless", "voiced")),
    conditions = list(c2place.ord = "velar")
)
## Summary:
##  * X : numeric predictor; with 100 values ranging from -34.964200 to 70.078500. 
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * c2place.ord : factor; set to the value(s): velar. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * rec.date : factor; set to the value(s): 24/03/2017 14:44:40.

3.1.3 PL04

filter(max, speaker == "pl04") %>%
    ggplot(aes(X, Y, colour = c2phonation)) +
    geom_point(alpha = 0.5) +
    facet_grid(vowel ~ c2place)

filter(max, speaker == "pl04", X > -30) %>%
    ggplot(aes(X, Y, colour = c2phonation)) +
    geom_point(alpha = 0.5) +
    facet_grid(vowel ~ c2place)

pl04_max <- filter(max, speaker == "pl04", X > -30)

pl04.m1 <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = pl04_max,
    method = "fREML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
rho <- start_value_rho(pl04.m1)

pl04.m1.ar <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = pl04_max,
    method = "ML",
    rho = rho,
    AR.start = pl04_max$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
pl04.m1.ar.null <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
#        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = pl04_max,
    method = "ML",
    rho = rho,
    AR.start = pl04_max$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
summary(pl04.m1.ar)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ X + s(X, bs = "cr") + s(X, by = c2phonation.ord, bs = "cr") + 
##     s(X, by = c2place.ord, bs = "cr") + s(X, by = vowel.ord, 
##     bs = "cr") + s(X, rec.date, bs = "fs", xt = "cr", m = 1, 
##     k = 5)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.28296    0.38386   8.553   <2e-16 ***
## X           -0.03567    0.03469  -1.028    0.304    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                edf  Ref.df      F  p-value    
## s(X)                         7.538   7.847 31.856  < 2e-16 ***
## s(X):c2phonation.ordvoiced   1.006   1.009  0.296    0.589    
## s(X):c2place.ordvelar        8.539   8.876 26.159  < 2e-16 ***
## s(X):vowel.ordo              7.975   8.620 10.694 7.37e-15 ***
## s(X,rec.date)              197.731 230.000 25.617  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Rank: 272/273
## R-sq.(adj) =  0.975   Deviance explained =   98%
## -ML = 1757.9  Scale est. = 0.72325   n = 1228
pl04.m1.null <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = pl04_max,
    method = "fREML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
# compareML(pl04.m1, pl04.m1.null)
compareML(pl04.m1.ar, pl04.m1.ar.null)
## pl04.m1.ar: Y ~ X + s(X, bs = "cr") + s(X, by = c2phonation.ord, bs = "cr") + 
##     s(X, by = c2place.ord, bs = "cr") + s(X, by = vowel.ord, 
##     bs = "cr") + s(X, rec.date, bs = "fs", xt = "cr", m = 1, 
##     k = 5)
## 
## pl04.m1.ar.null: Y ~ X + s(X, bs = "cr") + s(X, by = c2place.ord, bs = "cr") + 
##     s(X, by = vowel.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5)
## 
## Chi-square test of ML scores
## -----
##             Model    Score Edf Chisq    Df p.value Sig.
## 1 pl04.m1.ar.null 1758.008  11                         
## 2      pl04.m1.ar 1757.861  13 0.148 2.000   0.863
## Warning in compareML(pl04.m1.ar, pl04.m1.ar.null): AIC is not reliable,
## because an AR1 model is included (rho1 = 0.518140, rho2 = 0.518140).
## Warning in compareML(pl04.m1.ar, pl04.m1.ar.null): Only small difference in ML...
plot_gamsd(
    pl04.m1.ar,
    view = "X",
    comparison = list(c2phonation.ord = c("voiceless", "voiced")),
    conditions = list(c2place.ord = "coronal")
)
## Summary:
##  * X : numeric predictor; with 100 values ranging from -29.959700 to 32.556800. 
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * rec.date : factor; set to the value(s): 26/05/2017 17:58:11.

ggsave("presentations/2017-ultrafest/pl04.pdf", width = 6, height = 4)
acf_plot(resid(pl04.m1.ar), split_by=list(pl04_max$rec.date))

acf_resid(pl04.m1.ar, split_pred = "AR.start")

pl04_a <- filter(max, speaker == "pl04", vowel == "a")

pl04.m1 <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = pl04_a,
    method = "fREML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
rho <- start_value_rho(pl04.m1)

pl04.m1.ar <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = pl04_a,
    method = "ML",
    rho = rho,
    AR.start = pl04_a$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
pl04.m1.ar.null <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
#        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = pl04_a,
    method = "ML",
    rho = rho,
    AR.start = pl04_a$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
summary(pl04.m1.ar)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ X + s(X, bs = "cr") + s(X, by = c2phonation.ord, bs = "cr") + 
##     s(X, by = c2place.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  0.83223    0.46029   1.808   0.0711 .
## X            0.02144    0.04680   0.458   0.6470  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                               edf  Ref.df      F p-value    
## s(X)                        7.479   7.786 24.294  <2e-16 ***
## s(X):c2phonation.ordvoiced  3.384   4.190  0.553   0.575    
## s(X):c2place.ordvelar       8.019   8.623 12.720  <2e-16 ***
## s(X,rec.date)              95.926 116.000 29.624  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Rank: 148/149
## R-sq.(adj) =  0.987   Deviance explained = 98.9%
## -ML =   1101  Scale est. = 0.95795   n = 689
compareML(pl04.m1.ar, pl04.m1.ar.null)
## pl04.m1.ar: Y ~ X + s(X, bs = "cr") + s(X, by = c2phonation.ord, bs = "cr") + 
##     s(X, by = c2place.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5)
## 
## pl04.m1.ar.null: Y ~ X + s(X, bs = "cr") + s(X, by = c2place.ord, bs = "cr") + 
##     s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5)
## 
## Chi-square test of ML scores
## -----
##             Model    Score Edf Chisq    Df p.value Sig.
## 1 pl04.m1.ar.null 1101.707   9                         
## 2      pl04.m1.ar 1101.011  11 0.696 2.000   0.499
## Warning in compareML(pl04.m1.ar, pl04.m1.ar.null): AIC is not reliable,
## because an AR1 model is included (rho1 = 0.455701, rho2 = 0.455701).
## Warning in compareML(pl04.m1.ar, pl04.m1.ar.null): Only small difference in ML...

If only /a/ is used, no significance at all. With /o/ and X > -30, no significance, but with all X some significant difference at the very left end of the tongue.

plot_gamsd(
    pl04.m1.ar,
    view = "X",
    comparison = list(c2phonation.ord = c("voiceless", "voiced")),
    conditions = list(c2place.ord = "coronal")
)
## Summary:
##  * X : numeric predictor; with 100 values ranging from -41.892500 to 32.556800. 
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * rec.date : factor; set to the value(s): 26/05/2017 17:59:49.

3.1.4 PL03-PL04

pl34.m1 <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5) +
        s(X, speaker, bs = "fs", xt = "cr", m = 1, k = 5),
    data = max_pl_34,
    method = "fREML"
)

rho <- start_value_rho(pl34.m1)

pl34.m1.ar <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5) +
        s(X, speaker, bs = "fs", xt = "cr", m = 1, k = 5),
    data = max_pl_34,
    method = "ML",
    rho = rho,
    AR.start = max_pl_34$start.event
)

pl34.m1.ar.null <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
#        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5) +
        s(X, speaker, bs = "fs", xt = "cr", m = 1, k = 5),
    data = max_pl_34,
    method = "ML",
    rho = rho,
    AR.start = max_pl_34$start.event
)

summary(pl34.m1.ar)
compareML(pl34.m1.ar, pl34.m1.ar.null)
plot_gamsd(
    pl34.m1.ar,
    view = "X",
    comparison = list(c2phonation.ord = c("voiceless", "voiced")),
    conditions = list(c2place.ord = "velar", speaker = "pl04")
)

The models with X, Y are not reliable because the tongue size differs a lot between speakers.

pl34.m1.re <- bam(
    Y.re ~
        X.re +
        s(X.re, bs = "cr") +
        s(X.re, by = c2phonation.ord, bs = "cr") +
        s(X.re, by = c2place.ord, bs = "cr") +
        s(X.re, by = vowel.ord, bs = "cr") +
        s(X.re, rec.date, bs = "fs", xt = "cr", m = 1, k = 5) +
        s(X.re, speaker, bs = "fs", xt = "cr", m = 1, k = 5),
    data = max_pl_34,
    method = "fREML"
)

rho <- start_value_rho(pl34.m1.re)

pl34.m1.ar.re <- bam(
    Y.re ~
        X.re +
        s(X.re, bs = "cr") +
        s(X.re, by = c2phonation.ord, bs = "cr") +
        s(X.re, by = c2place.ord, bs = "cr") +
        s(X.re, by = vowel.ord, bs = "cr") +
        s(X.re, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = max_pl_34,
    method = "ML",
    rho = rho,
    AR.start = max_pl_34$start.event
)

pl34.m1.ar.re.null <- bam(
    Y.re ~
        X.re +
        s(X.re, bs = "cr") +
#        s(X, by = c2phonation.ord, bs = "cr") +
        s(X.re, by = c2place.ord, bs = "cr") +
        s(X.re, by = vowel.ord, bs = "cr") +
        s(X.re, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = max_pl_34,
    method = "ML",
    rho = rho,
    AR.start = max_pl_34$start.event
)

summary(pl34.m1.ar.re)
compareML(pl34.m1.ar.re, pl34.m1.ar.re.null)
plot_gamsd(
    pl34.m1.ar.re,
    view = "X.re",
    comparison = list(c2phonation.ord = c("voiceless", "voiced")),
    conditions = list(c2place.ord = "coronal")
)
plot_gamsd(
    pl34.m1.ar.re,
    view = "X.re",
    comparison = list(c2phonation.ord = c("voiceless", "voiced")),
    conditions = list(c2place.ord = "velar")
)

3.2 Italian

it01 and it02 have TRA at maximum displacement.

it.m1 <- bam(
    Y.re ~
        X.re +
        s(X.re, bs = "cr") +
        s(X.re, by = c2phonation.ord, bs = "cr") +
        s(X.re, by = c2place.ord, bs = "cr") +
        s(X.re, by = vowel.ord, bs = "cr") +
        s(X.re, rec.date, bs = "fs", xt = "cr", m = 1, k = 5) +
        s(X.re, speaker, bs = "fs", xt = "cr", m = 1, k = 5),
    data = max_it_12,
    method = "fREML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
rho <- start_value_rho(it.m1)

it.m1.ar <- bam(
    Y.re ~
        X.re +
        s(X.re, bs = "cr") +
        s(X.re, by = c2phonation.ord, bs = "cr") +
        s(X.re, by = c2place.ord, bs = "cr") +
        s(X.re, by = vowel.ord, bs = "cr") +
        s(X.re, rec.date, bs = "fs", xt = "cr", m = 1, k = 5) +
        s(X.re, speaker, bs = "fs", xt = "cr", m = 1, k = 5),
    data = max_it_12,
    method = "ML",
    rho = rho,
    AR.start = max_it_12$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
summary(it.m1.ar)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y.re ~ X.re + s(X.re, bs = "cr") + s(X.re, by = c2phonation.ord, 
##     bs = "cr") + s(X.re, by = c2place.ord, bs = "cr") + s(X.re, 
##     by = vowel.ord, bs = "cr") + s(X.re, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5) + s(X.re, speaker, bs = "fs", xt = "cr", 
##     m = 1, k = 5)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.25090    0.06259   4.008 6.26e-05 ***
## X.re         0.80000    0.21532   3.715 0.000207 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                   edf  Ref.df       F  p-value    
## s(X.re)                         5.516   6.082   3.409  0.00394 ** 
## s(X.re):c2phonation.ordvoiced   5.106   6.209   5.491 1.33e-05 ***
## s(X.re):c2place.ordvelar        8.706   8.942 101.743  < 2e-16 ***
## s(X.re):vowel.ordo              5.286   6.405   8.554 1.45e-09 ***
## s(X.re,rec.date)              177.731 415.000   1.097  < 2e-16 ***
## s(X.re,speaker)                 6.607   8.000  74.894  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Rank: 467/468
## R-sq.(adj) =  0.941   Deviance explained = 94.5%
## -ML = -6713.9  Scale est. = 0.0017165  n = 3165
it.m1.ar.null <- bam(
    Y.re ~
        X.re +
        s(X.re, bs = "cr") +
#        s(X.re, by = c2phonation.ord, bs = "cr") +
        s(X.re, by = c2place.ord, bs = "cr") +
        s(X.re, by = vowel.ord, bs = "cr") +
        s(X.re, rec.date, bs = "fs", xt = "cr", m = 1, k = 5) +
        s(X.re, speaker, bs = "fs", xt = "cr", m = 1, k = 5),
    data = max_it_12,
    method = "ML",
    rho = rho,
    AR.start = max_it_12$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
it.m1.null <- bam(
    Y.re ~
        X.re +
        s(X.re, bs = "cr") +
        s(X.re, by = c2place.ord, bs = "cr") +
        s(X.re, by = vowel.ord, bs = "cr") +
        s(X.re, rec.date, bs = "fs", xt = "cr", m = 1, k = 5) +
        s(X.re, speaker, bs = "fs", xt = "cr", m = 1, k = 5),
    data = max_it_12,
    method = "ML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
compareML(it.m1.ar, it.m1.ar.null)
## it.m1.ar: Y.re ~ X.re + s(X.re, bs = "cr") + s(X.re, by = c2phonation.ord, 
##     bs = "cr") + s(X.re, by = c2place.ord, bs = "cr") + s(X.re, 
##     by = vowel.ord, bs = "cr") + s(X.re, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5) + s(X.re, speaker, bs = "fs", xt = "cr", 
##     m = 1, k = 5)
## 
## it.m1.ar.null: Y.re ~ X.re + s(X.re, bs = "cr") + s(X.re, by = c2place.ord, 
##     bs = "cr") + s(X.re, by = vowel.ord, bs = "cr") + s(X.re, 
##     rec.date, bs = "fs", xt = "cr", m = 1, k = 5) + s(X.re, speaker, 
##     bs = "fs", xt = "cr", m = 1, k = 5)
## 
## Chi-square test of ML scores
## -----
##           Model     Score Edf  Chisq    Df   p.value Sig.
## 1 it.m1.ar.null -6702.103  14                            
## 2      it.m1.ar -6713.863  16 11.761 2.000 7.805e-06  ***
## Warning in compareML(it.m1.ar, it.m1.ar.null): AIC is not reliable, because
## an AR1 model is included (rho1 = 0.762798, rho2 = 0.762798).
plot_gamsd(
    it.m1.ar,
    view = "X.re",
    comparison = list(c2phonation.ord = c("voiceless", "voiced")),
    conditions = list(c2place.ord = "coronal")
)
## Summary:
##  * X.re : numeric predictor; with 100 values ranging from 0.000000 to 0.995850. 
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * rec.date : factor; set to the value(s): 29/11/2016 15:10:52. 
##  * speaker : factor; set to the value(s): it01.

3.2.1 IT01

it01_max <- filter(max, speaker == "it01")

it01.m1 <- bam(
    Y ~
        X.re +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = it01_max,
    method = "fREML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
rho <- start_value_rho(it01.m1)

it01.m1.ar <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = it01_max,
    method = "ML",
    rho = rho,
    AR.start = it01_max$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
summary(it01.m1.ar)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ X + s(X, bs = "cr") + s(X, by = c2phonation.ord, bs = "cr") + 
##     s(X, by = c2place.ord, bs = "cr") + s(X, by = vowel.ord, 
##     bs = "cr") + s(X, rec.date, bs = "fs", xt = "cr", m = 1, 
##     k = 5)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -4.73141    0.30682  -15.42   <2e-16 ***
## X            0.76532    0.02666   28.70   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                               edf  Ref.df       F  p-value    
## s(X)                        7.764   7.944 181.554  < 2e-16 ***
## s(X):c2phonation.ordvoiced  5.037   6.136   9.480 2.06e-10 ***
## s(X):c2place.ordvelar       8.822   8.979 183.531  < 2e-16 ***
## s(X):vowel.ordo             6.862   7.808  11.566 1.19e-15 ***
## s(X,rec.date)              92.680 225.000   1.976  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Rank: 267/268
## R-sq.(adj) =  0.979   Deviance explained =   98%
## -ML = 3488.2  Scale est. = 4.0701    n = 1932
it01.m1.ar.null <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
#        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = it01_max,
    method = "ML",
    rho = rho,
    AR.start = it01_max$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
compareML(it01.m1.ar, it01.m1.ar.null)
## it01.m1.ar: Y ~ X + s(X, bs = "cr") + s(X, by = c2phonation.ord, bs = "cr") + 
##     s(X, by = c2place.ord, bs = "cr") + s(X, by = vowel.ord, 
##     bs = "cr") + s(X, rec.date, bs = "fs", xt = "cr", m = 1, 
##     k = 5)
## 
## it01.m1.ar.null: Y ~ X + s(X, bs = "cr") + s(X, by = c2place.ord, bs = "cr") + 
##     s(X, by = vowel.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5)
## 
## Chi-square test of ML scores
## -----
##             Model    Score Edf  Chisq    Df   p.value Sig.
## 1 it01.m1.ar.null 3509.547  11                            
## 2      it01.m1.ar 3488.231  13 21.316 2.000 5.527e-10  ***
## Warning in compareML(it01.m1.ar, it01.m1.ar.null): AIC is not reliable,
## because an AR1 model is included (rho1 = 0.745137, rho2 = 0.745137).
plot_gamsd(
    it01.m1.ar,
    view = "X",
    comparison = list(c2phonation.ord = c("voiceless", "voiced")),
    conditions = list(c2place.ord = "coronal")
)
## Summary:
##  * X : numeric predictor; with 100 values ranging from -40.731200 to 54.425500. 
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * rec.date : factor; set to the value(s): 29/11/2016 15:10:52.

ggsave("presentations/2017-ultrafest/it01.pdf", width = 6, height = 4)

3.2.2 IT02

it02_max <- filter(max, speaker == "it02")

it02.m1 <- bam(
    Y ~
        X.re +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = it02_max,
    method = "fREML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
rho <- start_value_rho(it02.m1)

it02.m1.ar <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = it02_max,
    method = "ML",
    rho = rho,
    AR.start = it02_max$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
summary(it02.m1.ar)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ X + s(X, bs = "cr") + s(X, by = c2phonation.ord, bs = "cr") + 
##     s(X, by = c2place.ord, bs = "cr") + s(X, by = vowel.ord, 
##     bs = "cr") + s(X, rec.date, bs = "fs", xt = "cr", m = 1, 
##     k = 5)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.67360    0.52437  -3.192  0.00145 ** 
## X            0.65159    0.03157  20.640  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                               edf  Ref.df      F  p-value    
## s(X)                        6.386   7.189 57.403  < 2e-16 ***
## s(X):c2phonation.ordvoiced  6.589   7.613  4.927 9.93e-06 ***
## s(X):c2place.ordvelar       8.630   8.910 52.382  < 2e-16 ***
## s(X):vowel.ordo             5.937   7.020  9.062 6.16e-11 ***
## s(X,rec.date)              72.997 185.000  1.401  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Rank: 227/228
## R-sq.(adj) =  0.916   Deviance explained = 92.2%
## -ML = 2803.6  Scale est. = 9.778     n = 1233
it02.m1.ar.null <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
#        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = it02_max,
    method = "ML",
    rho = rho,
    AR.start = it02_max$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
compareML(it02.m1.ar, it02.m1.ar.null)
## it02.m1.ar: Y ~ X + s(X, bs = "cr") + s(X, by = c2phonation.ord, bs = "cr") + 
##     s(X, by = c2place.ord, bs = "cr") + s(X, by = vowel.ord, 
##     bs = "cr") + s(X, rec.date, bs = "fs", xt = "cr", m = 1, 
##     k = 5)
## 
## it02.m1.ar.null: Y ~ X + s(X, bs = "cr") + s(X, by = c2place.ord, bs = "cr") + 
##     s(X, by = vowel.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5)
## 
## Chi-square test of ML scores
## -----
##             Model    Score Edf Chisq    Df   p.value Sig.
## 1 it02.m1.ar.null 2813.488  11                           
## 2      it02.m1.ar 2803.592  13 9.897 2.000 5.034e-05  ***
## Warning in compareML(it02.m1.ar, it02.m1.ar.null): AIC is not reliable,
## because an AR1 model is included (rho1 = 0.735050, rho2 = 0.735050).
plot_gamsd(
    it02.m1.ar,
    view = "X",
    comparison = list(c2phonation.ord = c("voiceless", "voiced")),
    conditions = list(c2place.ord = "coronal")
)
## Summary:
##  * X : numeric predictor; with 100 values ranging from -44.949000 to 62.666500. 
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * rec.date : factor; set to the value(s): 12/12/2016 14:45:14.

4 Italian tongue at closure

tongues_clos %>%
    filter(c2place == "coronal", language == "italian") %>%
    plot_tongue(geom = "point", alpha = 0.5) +
    aes(colour = c2phonation) +
    facet_grid(vowel ~ speaker)

4.1 IT01

it01_clos <- tongues_clos %>%
    filter(speaker == "it01", vowel != "u")
it01_clos_m1 <- bam(
    Y ~
        X.re +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = it01_clos,
    method = "fREML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
rho <- start_value_rho(it01_clos_m1)

it01_clos_m1_ar <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = it01_clos,
    method = "ML",
    rho = rho,
    AR.start = it01_clos$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
summary(it01_clos_m1_ar)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ X + s(X, bs = "cr") + s(X, by = c2phonation.ord, bs = "cr") + 
##     s(X, by = c2place.ord, bs = "cr") + s(X, by = vowel.ord, 
##     bs = "cr") + s(X, rec.date, bs = "fs", xt = "cr", m = 1, 
##     k = 5)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.61066    0.39198   1.558    0.119    
## X            0.84888    0.02509  33.831   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                               edf  Ref.df       F  p-value    
## s(X)                        7.549   7.867 198.722  < 2e-16 ***
## s(X):c2phonation.ordvoiced  4.683   5.594  13.824 3.56e-14 ***
## s(X):c2place.ordvelar       8.641   8.924 174.974  < 2e-16 ***
## s(X):vowel.ordo             6.049   6.977  14.770  < 2e-16 ***
## s(X,rec.date)              96.251 230.000   1.789  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Rank: 272/273
## R-sq.(adj) =  0.981   Deviance explained = 98.2%
## -ML = 2983.9  Scale est. = 3.7241    n = 1648
it01_clos_m1_ar_null <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
#        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = it01_clos,
    method = "ML",
    rho = rho,
    AR.start = it01_clos$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
compareML(it01_clos_m1_ar, it01_clos_m1_ar_null)
## it01_clos_m1_ar: Y ~ X + s(X, bs = "cr") + s(X, by = c2phonation.ord, bs = "cr") + 
##     s(X, by = c2place.ord, bs = "cr") + s(X, by = vowel.ord, 
##     bs = "cr") + s(X, rec.date, bs = "fs", xt = "cr", m = 1, 
##     k = 5)
## 
## it01_clos_m1_ar_null: Y ~ X + s(X, bs = "cr") + s(X, by = c2place.ord, bs = "cr") + 
##     s(X, by = vowel.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5)
## 
## Chi-square test of ML scores
## -----
##                  Model    Score Edf  Chisq    Df   p.value Sig.
## 1 it01_clos_m1_ar_null 3014.072  11                            
## 2      it01_clos_m1_ar 2983.854  13 30.218 2.000 7.523e-14  ***
## Warning in compareML(it01_clos_m1_ar, it01_clos_m1_ar_null): AIC is
## not reliable, because an AR1 model is included (rho1 = 0.721734, rho2 =
## 0.721734).
plot_gamsd(
    it01_clos_m1_ar,
    view = "X",
    comparison = list(c2phonation.ord = c("voiceless", "voiced")),
    conditions = list(c2place.ord = "coronal")
)
## Summary:
##  * X : numeric predictor; with 100 values ranging from -39.125800 to 54.502000. 
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * rec.date : factor; set to the value(s): 29/11/2016 15:10:52.

ggsave("presentations/2017-ultrafest/it01-clos.pdf", width = 6, height = 4)

4.2 IT02

it02_clos <- tongues_clos %>%
    filter(speaker == "it02", vowel != "u")
it02_clos_m1 <- bam(
    Y ~
        X.re +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = it02_clos,
    method = "fREML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
rho <- start_value_rho(it02_clos_m1)

it02_clos_m1_ar <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = it02_clos,
    method = "ML",
    rho = rho,
    AR.start = it02_clos$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
summary(it02_clos_m1_ar)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ X + s(X, bs = "cr") + s(X, by = c2phonation.ord, bs = "cr") + 
##     s(X, by = c2place.ord, bs = "cr") + s(X, by = vowel.ord, 
##     bs = "cr") + s(X, rec.date, bs = "fs", xt = "cr", m = 1, 
##     k = 5)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.56479    0.48605  -5.277 1.57e-07 ***
## X            0.64094    0.04079  15.714  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                edf  Ref.df      F  p-value    
## s(X)                         7.047   7.612 40.215  < 2e-16 ***
## s(X):c2phonation.ordvoiced   7.191   8.131  5.061 2.41e-06 ***
## s(X):c2place.ordvelar        8.088   8.702 29.725  < 2e-16 ***
## s(X):vowel.ordo              8.257   8.776 10.910 2.49e-15 ***
## s(X,rec.date)              115.667 195.000  2.711  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Rank: 237/238
## R-sq.(adj) =   0.92   Deviance explained = 92.8%
## -ML = 2848.9  Scale est. = 6.7751    n = 1319
it02_clos_m1_ar_null <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
#        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = it02_clos,
    method = "ML",
    rho = rho,
    AR.start = it02_clos$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
compareML(it02_clos_m1_ar, it02_clos_m1_ar_null)
## it02_clos_m1_ar: Y ~ X + s(X, bs = "cr") + s(X, by = c2phonation.ord, bs = "cr") + 
##     s(X, by = c2place.ord, bs = "cr") + s(X, by = vowel.ord, 
##     bs = "cr") + s(X, rec.date, bs = "fs", xt = "cr", m = 1, 
##     k = 5)
## 
## it02_clos_m1_ar_null: Y ~ X + s(X, bs = "cr") + s(X, by = c2place.ord, bs = "cr") + 
##     s(X, by = vowel.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5)
## 
## Chi-square test of ML scores
## -----
##                  Model    Score Edf  Chisq    Df   p.value Sig.
## 1 it02_clos_m1_ar_null 2861.022  11                            
## 2      it02_clos_m1_ar 2848.944  13 12.078 2.000 5.680e-06  ***
## Warning in compareML(it02_clos_m1_ar, it02_clos_m1_ar_null): AIC is
## not reliable, because an AR1 model is included (rho1 = 0.718876, rho2 =
## 0.718876).
plot_gamsd(
    it02_clos_m1_ar,
    view = "X",
    comparison = list(c2phonation.ord = c("voiceless", "voiced")),
    conditions = list(c2place.ord = "coronal")
)
## Summary:
##  * X : numeric predictor; with 100 values ranging from -45.980000 to 64.280100. 
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): o. 
##  * rec.date : factor; set to the value(s): 12/12/2016 14:44:52.

Both it01 and it02 show TRA at closure. It is difficult to quantify the difference compared to TRA at maximum displacement. Need a model to check that.

5 Polish tongue at closure

5.1 PL04

pl04_clos <- tongues_clos %>%
    filter(speaker == "pl04", vowel != "u")
pl04_clos_m1 <- bam(
    Y ~
        X.re +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = pl04_clos,
    method = "fREML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
rho <- start_value_rho(pl04_clos_m1)

pl04_clos_m1_ar <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = pl04_clos,
    method = "ML",
    rho = rho,
    AR.start = pl04_clos$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
summary(pl04_clos_m1_ar)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Y ~ X + s(X, bs = "cr") + s(X, by = c2phonation.ord, bs = "cr") + 
##     s(X, by = c2place.ord, bs = "cr") + s(X, by = vowel.ord, 
##     bs = "cr") + s(X, rec.date, bs = "fs", xt = "cr", m = 1, 
##     k = 5)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.53763    0.11242 -13.677   <2e-16 ***
## X           -0.01176    0.01746  -0.674    0.501    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                edf  Ref.df       F  p-value    
## s(X)                         7.817   7.951 234.521  < 2e-16 ***
## s(X):c2phonation.ordvoiced   1.004   1.005   0.809    0.369    
## s(X):c2place.ordvelar        5.960   7.099   5.674 1.66e-06 ***
## s(X):vowel.ordo              7.493   8.406  16.701  < 2e-16 ***
## s(X,rec.date)              155.677 234.000   5.321  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Rank: 277/278
## R-sq.(adj) =  0.984   Deviance explained = 98.6%
## -ML = 1886.4  Scale est. = 0.83789   n = 1387
pl04_clos_m1_ar_null <- bam(
    Y ~
        X +
        s(X, bs = "cr") +
#        s(X, by = c2phonation.ord, bs = "cr") +
        s(X, by = c2place.ord, bs = "cr") +
        s(X, by = vowel.ord, bs = "cr") +
        s(X, rec.date, bs = "fs", xt = "cr", m = 1, k = 5),
    data = pl04_clos,
    method = "ML",
    rho = rho,
    AR.start = pl04_clos$start.event
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
compareML(pl04_clos_m1_ar, pl04_clos_m1_ar_null)
## pl04_clos_m1_ar: Y ~ X + s(X, bs = "cr") + s(X, by = c2phonation.ord, bs = "cr") + 
##     s(X, by = c2place.ord, bs = "cr") + s(X, by = vowel.ord, 
##     bs = "cr") + s(X, rec.date, bs = "fs", xt = "cr", m = 1, 
##     k = 5)
## 
## pl04_clos_m1_ar_null: Y ~ X + s(X, bs = "cr") + s(X, by = c2place.ord, bs = "cr") + 
##     s(X, by = vowel.ord, bs = "cr") + s(X, rec.date, bs = "fs", 
##     xt = "cr", m = 1, k = 5)
## 
## Chi-square test of ML scores
## -----
##                  Model    Score Edf Chisq    Df p.value Sig.
## 1 pl04_clos_m1_ar_null 1886.770  11                         
## 2      pl04_clos_m1_ar 1886.368  13 0.401 2.000   0.669
## Warning in compareML(pl04_clos_m1_ar, pl04_clos_m1_ar_null): AIC is
## not reliable, because an AR1 model is included (rho1 = 0.494196, rho2 =
## 0.494196).
## Warning in compareML(pl04_clos_m1_ar, pl04_clos_m1_ar_null): Only small difference in ML...
plot_gamsd(
    pl04_clos_m1_ar,
    view = "X",
    comparison = list(c2phonation.ord = c("voiceless", "voiced")),
    conditions = list(c2place.ord = "coronal")
)
## Summary:
##  * X : numeric predictor; with 100 values ranging from -39.613600 to 30.560600. 
##  * c2phonation.ord : factor; set to the value(s): voiced, voiceless. 
##  * c2place.ord : factor; set to the value(s): coronal. 
##  * vowel.ord : factor; set to the value(s): a. 
##  * rec.date : factor; set to the value(s): 26/05/2017 17:59:34.

ggsave("presentations/2017-ultrafest/pl04-clos.pdf", width = 6, height = 4)